{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "e6b6aa0f",
   "metadata": {},
   "source": [
    "# Newton-Cotes\n",
    "1. 本质上是使用积分中值定理+线性插值公式\n",
    "2. 根据选取插值点的数量分为不同类型的数值积分公式\n",
    "3. <font color=maroon><b>newton_cotes()</b></font>函数将返回差值点全值以及相应的余项"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "4f4f5555",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([ 0.27908289,  1.66151675, -0.26186949,  2.96183422, -1.28112875,\n",
       "         2.96183422, -0.26186949,  1.66151675,  0.27908289]),\n",
       " -0.0050622628400406175)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy import integrate\n",
    "    \n",
    "integrate.newton_cotes(8,equal=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fa0229da",
   "metadata": {},
   "source": [
    "# Romberg\n",
    "1. 是复化积分公式的加速版本\n",
    "2. 相比较于独立的复化N-C公式而言，取更少的分点，计算时有更少的四则运算\n",
    "3. <font color=maroon><b>romberg(func,a,b,args=(),tol,show=False)</b></font>能够返回积分计算结果\n",
    "    > 1. tol是最小误差<br>\n",
    "    > 2. show打印推导过程<br>\n",
    "4. <font color=maroon><b>romb(y,dx,axis=-1,show=False)</b></font>同样能够返回积分计算结果，但是这个函数完成的功能是直接根据给出的函数值作为$f(x_k)$进行计算\n",
    "    > 1. axis是积分关注的维度<br>\n",
    "    > 2. dx调整积分插值点之间的步长，一般情况下有$(2^k+1)$个数据点的情况下步长应为$\\frac{(b-a)}{k}$<br>\n",
    "    > 2. show打印推导过程<br>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6e0d9a1d",
   "metadata": {},
   "source": [
    "$$\n",
    "    \\displaystyle\\int_{0}^{1}{\\frac{4}{1+x^2}dx}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "e5628e67",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Romberg integration of <function vectorize1.<locals>.vfunc at 0x7fbc2bc7cc10> from [0, 1]\n",
      "\n",
      " Steps  StepSize   Results\n",
      "     1  1.000000  3.000000 \n",
      "     2  0.500000  3.100000  3.133333 \n",
      "     4  0.250000  3.131176  3.141569  3.142118 \n",
      "     8  0.125000  3.138988  3.141593  3.141594  3.141586 \n",
      "    16  0.062500  3.140942  3.141593  3.141593  3.141593  3.141593 \n",
      "    32  0.031250  3.141430  3.141593  3.141593  3.141593  3.141593  3.141593 \n",
      "\n",
      "The final result is 3.141592653638244 after 33 function evaluations.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "3.141592653638244"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy import integrate\n",
    "\n",
    "def func(x):\n",
    "    return 4/(1 + x**2)\n",
    "resx = integrate.romberg(func, 0, 1, show=True)\n",
    "resx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "baf4e30d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "       Richardson Extrapolation Table for Romberg Integration       \n",
      "====================================================================\n",
      " 3.00000 \n",
      " 3.10000  3.13333 \n",
      " 3.13118  3.14157  3.14212 \n",
      " 3.13899  3.14159  3.14159  3.14159 \n",
      " 3.14094  3.14159  3.14159  3.14159  3.14159 \n",
      " 3.14143  3.14159  3.14159  3.14159  3.14159  3.14159 \n",
      " 3.14155  3.14159  3.14159  3.14159  3.14159  3.14159  3.14159 \n",
      "====================================================================\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "3.141592653589722"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy import integrate\n",
    "import numpy as np\n",
    "\n",
    "x = np.linspace(0,1,65,endpoint=True)\n",
    "y = 4/(1 + np.power(x,2))\n",
    "resx = integrate.romb(y,dx=1/64,show=True)\n",
    "resx"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ba36c51c",
   "metadata": {},
   "source": [
    "# Gaussian\n",
    "1. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1dbae710",
   "metadata": {},
   "source": [
    "# python implementation\n",
    "## single integral\n",
    "1. quad(func,a,b)\n",
    "2. quad_vec(func,a,b)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4dd7a801",
   "metadata": {},
   "source": [
    "### sample1\n",
    "$$\n",
    "    {x}^2 + {y}^2 = 1\n",
    "$$\n",
    "$$\n",
    "    y = \\sqrt{1-x^2}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "10994797",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "res of scipy quad\n",
      "1.5707963267948983 1.0002354500215915e-09\n",
      "res of closed formula\n",
      "1.5707963267948966 1.7763568394002505e-15\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "from scipy import constants\n",
    "from scipy import integrate\n",
    "\n",
    "# calculate area of half circle\n",
    "def half_circle(x):\n",
    "    return (1-x**2)**0.5\n",
    "\n",
    "resx, err = integrate.quad(half_circle, -1, 1)\n",
    "print(\"res of scipy quad\")\n",
    "print(resx, err)\n",
    "print(\"res of closed formula\")\n",
    "print(constants.pi *0.5, abs(constants.pi *0.5 - resx))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7437eef7",
   "metadata": {},
   "source": [
    "## double integral & triple integral\n",
    "1. dblquad(func2d, a,b, gfun, hfun, args=())\n",
    "2. tblquad(func2d, a,b, gfun, hfun, qfun ,rfun, args=())\n",
    "\n",
    "其中，$(a,b)$是x的积分上下限，$(gfun,hfun)$是$y(x)$的积分上下限。args是传入func2d中的变量数值"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "732ea5ae",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "python3.6",
   "language": "python",
   "name": "python3.6"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {
    "height": "calc(100% - 180px)",
    "left": "10px",
    "top": "150px",
    "width": "315px"
   },
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
